Objectif : Exploration des données et modélisations
Date début de l’analyse : 17 mai 2022
Date de la dernière modification : 24 mai 2022
# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts", "DescTools")
## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)}})
# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))
# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
oa_color_openalex == "gold" ~ "gold",
oa_color_openalex == "hybrid" ~ "hybridOA",
TRUE ~ NA_character_), #qs Maurits : catégorie "other" ou tout en NA ? Sur quelle variable se baser alors ?
amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
TRUE ~ amount_apc_EUR))
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup()
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%
group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>%
mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :
Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.
# Sélection de variables pour les modèles
bso_mvp <- data %>%
select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>%
mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))
# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>%
filter(!is.na(bso_classification),
!is.na(oa_color_finale),
!is.na(is_french_CA_bso_wos_openalex_single_lang),
!is.na(amount_apc_EUR_trusted),
!is.na(is_covered_by_couperin)) %>%
filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)
À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :
apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.
graph <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
geom_boxplot(fill = "#112446") +
coord_flip() +
labs(x = "Montant des APC fiables", y = " ", title = "Boxplot de Y") +
theme_ipsum() +
theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12))
#boxplot <- ggplotly(graph)
graph1 <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
geom_histogram(bins = 30L, fill = "#112446") +
labs(x = "Montant des APC fiables", y = "Nombre d'articles", title = "Histogramme de Y") +
theme_ipsum() +
theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12)) +
stat_function(fun = dnorm, args = list(mean = mean(bso_pop$amount_apc_EUR_trusted), sd = sd(bso_pop$amount_apc_EUR_trusted)), col = "white")
#histo <- ggplotly(graph1)
grid.arrange(graph, graph1, ncol = 2, nrow = 1)
L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros, nous filtrons donc la base sur les montants inférieurs à 7000 euros (le 3e APC le plus élevé est de 6848.67 euros) et regardons à nouveau les graphiques d’évolution des APC par année, selon les différentes variables qualitatives.
# Suppression des valeurs aberrantes
bso_pop <- bso_pop %>% filter(amount_apc_EUR_trusted < 7000) %>% mutate(is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang))
# Graphique / discipline
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, bso_classification) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / tier
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / french CA
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>%
rename(is_french_CA = is_french_CA_bso_wos_openalex_single_lang)
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Graphique / oa_color_finale
# table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
# plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>%
summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))
# Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
geom_line(size=1, alpha=0.9, linetype=1) +
geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
#scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
theme_classic() +
theme(legend.position = "right")
ggplotly(graph)
# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
geom_histogram(aes(y=..density..), color="black", fill = "steelblue", alpha = 0.2) +
geom_density( color="red", size = 1, adjust = 5) +
stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd)) +
xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
xlab("Montant des APC fiables") + ylab("Densité") +
theme_classic()
# Test statistique de normalité
out <- JarqueBeraTest(bso_pop$amount_apc_EUR_trusted, robust = FALSE, method = "chisq") # Y non normal (on rejette H0)
library(formattable)
df <- data.frame(value = unlist(out))
tdf <- as.data.frame(t(df))
formattable(tdf)
| statistic.X-squared | parameter.df | p.value | method | data.name | |
|---|---|---|---|---|---|
| value | 47076.3232661981 | 2 | 0 | Jarque Bera Test | bso_pop$amount_apc_EUR_trusted |
Le montant des APC fiables n’est pas normalement distribué.
# Fonction pour sauvegarder les résultats des tests
save_results <- function(output, num){
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
assign(glue("tdf_{num}"), tdf[,1:5], envir = .GlobalEnv)
}
# Test de corrélation (Spearman car distribution non normale)
output <- cor.test(bso_pop$amount_apc_EUR_trusted, bso_pop$year, method = "spearman") # corrélation significative
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"), #escape special characters
data.name = gsub("bso_pop\\$", "", data.name))
formattable(tdf)
| statistic.S | p.value | estimate.rho | null.value.rho | alternative | method | data.name | |
|---|---|---|---|---|---|---|---|
| value | 242612563868318 | 0 | 0.149413478963885 | 0 | two.sided | Spearman’s rank correlation rho | amount_apc_EUR_trusted and year |
# Tests de dépendance
# Année
save_results(chisq.test(bso_pop$year, bso_pop$bso_classification), "1")
save_results(chisq.test(bso_pop$year, bso_pop$tier), "2")
save_results(chisq.test(bso_pop$year, bso_pop$oa_color_finale), "3")
save_results(chisq.test(bso_pop$year, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "4")
save_results(chisq.test(bso_pop$year, bso_pop$is_covered_by_couperin), "5")
# Discipline
save_results(chisq.test(bso_pop$bso_classification, bso_pop$tier), "6")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$oa_color_finale), "7")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "8")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_covered_by_couperin), "9")
# Tier
save_results(chisq.test(bso_pop$tier, bso_pop$oa_color_finale), "10")
save_results(chisq.test(bso_pop$tier, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "11")
save_results(chisq.test(bso_pop$tier, bso_pop$is_covered_by_couperin), "12")
# Couleur OA
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "13")
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_covered_by_couperin), "14")
# French CA
save_results(chisq.test(bso_pop$is_french_CA_bso_wos_openalex_single_lang, bso_pop$is_covered_by_couperin), "15") #toutes dépendantes
# Présentation des résultats dans une table
output <- rbind(tdf_1,tdf_2,tdf_3,tdf_4,tdf_5,tdf_6,tdf_7,tdf_8,tdf_9,tdf_10,tdf_11,tdf_12,tdf_13,tdf_14,tdf_15) %>%
mutate(p.value = as.numeric(p.value),
Dependance = case_when(p.value <= 0.05 ~ "Oui",
TRUE ~ "Non"))
output <- output %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"), #escape special characters
data.name = gsub("bso_pop\\$", "", data.name)) #remove base name
rownames(output) <- NULL #remove rownames
kable(output, format = "html", caption = "Résultats des tests de dépendance entre les variables qualitatives") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| statistic.X-squared | parameter.df | p.value | method | data.name | Dependance |
|---|---|---|---|---|---|
| 1843.61062820764 | 63 | 0 | Pearson’s Chi-squared test | year and bso_classification | Oui |
| 10563.0744012642 | 21 | 0 | Pearson’s Chi-squared test | year and tier | Oui |
| 89.3806015204807 | 7 | 0 | Pearson’s Chi-squared test | year and oa_color_finale | Oui |
| 207.610054207291 | 7 | 0 | Pearson’s Chi-squared test | year and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 155.89053408379 | 7 | 0 | Pearson’s Chi-squared test | year and is_covered_by_couperin | Oui |
| 7904.05509513946 | 27 | 0 | Pearson’s Chi-squared test | bso_classification and tier | Oui |
| 4499.27337899457 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and oa_color_finale | Oui |
| 223.812265757185 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 13108.8454392833 | 9 | 0 | Pearson’s Chi-squared test | bso_classification and is_covered_by_couperin | Oui |
| 11884.0898801833 | 3 | 0 | Pearson’s Chi-squared test | tier and oa_color_finale | Oui |
| 228.259385498635 | 3 | 0 | Pearson’s Chi-squared test | tier and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 15863.5696643588 | 3 | 0 | Pearson’s Chi-squared test | tier and is_covered_by_couperin | Oui |
| 1870.74427535178 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | oa_color_finale and is_french_CA_bso_wos_openalex_single_lang | Oui |
| 49845.1115001406 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | oa_color_finale and is_covered_by_couperin | Oui |
| 1145.57365987729 | 1 | 0 | Pearson’s Chi-squared test with Yates’ continuity correction | is_french_CA_bso_wos_openalex_single_lang and is_covered_by_couperin | Oui |
# Vérifier :
#- homoscédasticité (homogénéité de la variance au sein des sous-pops)
#- indépendance (sous-pops doivent avoir un effet aléatoire sur Y)
#- autocorrélation (résidus et erreurs ne doivent pas être auto-correlés)
Pour les modélisations, nous cherchons donc à expliquer et prédire le montant des APC fiables grâce aux variables suivantes :
# Train test split
train_size <- floor(0.8 * nrow(bso_pop)) # 80% of the sample size
## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(bso_pop)), size = train_size)
train <- bso_pop[train_ind, ]
test <- bso_pop[-train_ind, ]
Ensuite, nous séparons notre base de données en 2 échantillons : un échantillon d’apprentissage (train) et un échantillon de test (test). Le premier contient 80% des données, soit 95691 articles et le second contient les 20% derniers, soit 23923 articles. Tous les modèles seront construits sur l’échantillon d’apprentissage, puis seront validés / testés sur l’échantillon test ; s’agissant de données inconnues (les modèles n’auront pas été entraînés dessus), les appliquer à cet échantillon permettra de tester leur robustesse et leur capacité prédictive.
# Reg linéaire
lm1 <- lm(amount_apc_EUR_trusted ~ year, data = train)
summary <- tidy(lm1)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -77428.73049 | 2396.139753 | -32.31395 | 0 |
| year | 39.27705 | 1.187785 | 33.06747 | 0 |
Chaque année, le montant des APC (fiables) augmente de 39 euros en moyenne.
# Reg linéaire
lm2 <- lm(amount_apc_EUR_trusted ~ year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang, data = train)
summary <- tidy(lm2)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -60191.26359 | 2256.790367 | -26.67118 | 0 |
| year | 30.85978 | 1.118668 | 27.58619 | 0 |
| tier2 | -311.08861 | 9.693970 | -32.09094 | 0 |
| tier3 | -440.85370 | 9.574222 | -46.04590 | 0 |
| tier4 | -446.89060 | 10.263922 | -43.53995 | 0 |
| oa_color_finalehybridOA | 801.69500 | 8.181967 | 97.98316 | 0 |
| is_covered_by_couperin1 | -96.96469 | 9.262517 | -10.46850 | 0 |
| is_french_CA_bso_wos_openalex_single_lang1 | -67.89926 | 4.874324 | -13.92998 | 0 |
La discipline est ici la source de variations du modèle, chaque discipline dispose de son ordonnée à l’origine.
# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>%
set_engine("lmer")
# une ordonnée à l'origine qui peut varier d'une discipline à l'autre
model1_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ (1 | bso_classification), data = train)
# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1495.2871 | 105.4812 | 14.17586 |
| ran_pars | bso_classification | sd__(Intercept) | 332.7945 | NA | NA |
| ran_pars | Residual | sd__Observation | 810.2233 | NA | NA |
# les variations des APC par discipline
tidy(model1_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1777 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1777, linetype = 2) +
coord_flip()
Par rapport au premier modèle, on ajoute l’année comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par discipline).
## Tentons un 2e modèle, dans lequel on suppose que le montant des APC augmente au cours du temps
# on commence par transformer la variable année
train <- train %>%
mutate(year = year - 2013)
model2_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ year + (1 | bso_classification), data = train) # l'effet de l'année sur l'APC est supposé le même quel que soit la discipline
# vue générale du modèle
summary <- tidy(model2_fit) # l'APC moyen augmente de 34,6€ par an
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1311.20846 | 103.378135 | 12.68361 |
| fixed | NA | year | 41.67462 | 1.174069 | 35.49587 |
| ran_pars | bso_classification | sd__(Intercept) | 325.72589 | NA | NA |
| ran_pars | Residual | sd__Observation | 804.94652 | NA | NA |
# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model2_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1565 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1565, linetype = 2) +
coord_flip()
Par rapport au deuxième modèle, on autorise la pente liée à l’année, à varier selon la discipline (more random effects).
## Dans un 3e modèle, on suppose que le montant des APC augmente au cours du temps mais potentiellement de manière différente d'une discipline à l'autre
model3_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ year + (1 + year | bso_classification), data = train) # l'effet de l'année sur l'APC est dit aléatoire (il peut varier d'une discipline à l'autre)
# vue générale du modèle
summary <- tidy(model3_fit) # l'APC moyen augmente de 29.3€ par an, mais il y a pas mal de variabilité
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1439.5454903 | 99.17308 | 14.515486 |
| fixed | NA | year | 13.7111858 | 10.77209 | 1.272844 |
| ran_pars | bso_classification | sd__(Intercept) | 309.8740532 | NA | NA |
| ran_pars | bso_classification | cor__(Intercept).year | -0.0450307 | NA | NA |
| ran_pars | bso_classification | sd__year | 32.6505586 | NA | NA |
| ran_pars | Residual | sd__Observation | 802.6447803 | NA | NA |
# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model3_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
filter(term %in% "(Intercept)") %>% # dans un premier temps on ne regarde que l'AC moyen par discipline
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
mutate(estimate = 1925 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1925, linetype = 2) +
coord_flip()
tidy(model3_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
filter(term %in% "year") %>% # maintenant on ne garde que les pentes
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
mutate(estimate = 29.3 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 29.3, linetype = 2) +
coord_flip()
# dans les humanités ou en médecine, le montant des APC diminue au fil du temps tandis qu'en sciences sociales et en physique, il augmente
tidy(model3_fit, effects = "ran_vals") %>%
select(-`std.error`) %>%
pivot_wider(names_from = term, values_from = estimate)
effect group level (Intercept) year
L’année est ici la source de variations du modèle, chaque année dispose de son ordonnée à l’origine.
# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>%
set_engine("lmer")
# une ordonnée à l'origine qui peut varier d'une année à l'autre
model1_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ 1 + (1 | year), data = train)
# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1771.9066 | 41.84011 | 42.34947 |
| ran_pars | year | sd__(Intercept) | 118.0731 | NA | NA |
| ran_pars | Residual | sd__Observation | 817.8890 | NA | NA |
# les variations des APC par année
tidy(model1_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1772 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1772, linetype = 2) +
coord_flip()
Par rapport au premier modèle, on ajoute la discipline comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par année).
## Tentons un 2e modèle, dans lequel on suppose que le montant des APC dépend de la discipline
model2_fit <- model1_spec %>%
fit(amount_apc_EUR_trusted ~ bso_classification + (1 | year), data = train) # l'effet de la discipline sur l'APC est supposé le même quelle que soit l'année
# vue générale du modèle
summary <- tidy(model2_fit)
kable(summary, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = F,
html_font = "sans-serif")
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 1849.35886 | 43.361568 | 42.649723 |
| fixed | NA | bso_classificationChemistry | -197.16272 | 15.207234 | -12.965061 |
| fixed | NA | bso_classificationComputer and information sciences | -553.73092 | 19.908469 | -27.813838 |
| fixed | NA | bso_classificationEarth, Ecology, Energy and applied biology | -232.68578 | 9.258577 | -25.131917 |
| fixed | NA | bso_classificationEngineering | -486.09551 | 27.113104 | -17.928434 |
| fixed | NA | bso_classificationHumanities | -882.32097 | 34.474890 | -25.593148 |
| fixed | NA | bso_classificationMathematics | -925.75223 | 34.972179 | -26.471106 |
| fixed | NA | bso_classificationMedical research | -13.18786 | 6.038334 | -2.184023 |
| fixed | NA | bso_classificationPhysical sciences, Astronomy | -227.31944 | 10.420008 | -21.815669 |
| fixed | NA | bso_classificationSocial sciences | -417.68790 | 34.692059 | -12.039871 |
| ran_pars | year | sd__(Intercept) | 122.01518 | NA | NA |
| ran_pars | Residual | sd__Observation | 804.11607 | NA | NA |
# les variations des APC par année - ici les valeurs sont calculées pour la discipline "Social"
tidy(model2_fit, effects = "ran_vals") %>% # on calcule les valeurs des effets aléatoires
mutate(estimate = 1852 + estimate, # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
ymax = estimate + 2 * `std.error`) %>%
arrange(estimate) %>%
mutate(level = as_factor(level)) %>%
ggplot(aes(x = level, y = estimate)) +
geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 1852, linetype = 2) +
coord_flip()
Par rapport au deuxième modèle, on autorise la pente liée à la discipline, à varier selon l’année (more random effects).
### Hypothèses
# Test d'homogénéité de la variance (H0)
bartlett.test(amount_apc_EUR_trusted ~ interaction(year, is_covered_by_couperin), data = bso_pop)
Bartlett test of homogeneity of variances
data: amount_apc_EUR_trusted by interaction(year, is_covered_by_couperin) Bartlett’s K-squared = 1959.9, df = 15, p-value < 2.2e-16
#-->use weights arguments from lme()
# Test d'autocorrélation des résidus
# Taille des sous-populations
table(bso_pop$is_covered_by_couperin) # échantillons déséquilibrés
0 1
102667 16947
table(bso_pop$year) # de 8184 à 26054 (peu d'obs pour 2013-2014 car suppression des NA)
2013 2014 2015 2016 2017 2018 2019 2020 8184 9537 11251 12978 14839 16943 19828 26054
# 16 clusters
#--> passer en bayésien
# Filtre sur les Y is_french_CA manquants
sample <- data %>% filter(is.na(is_french_CA_bso_wos_openalex_single_lang))
Document sous licence ouverte